Rank-normalized split-Rhat and relative efficiency estimates
Aki Vehtari, Daniel Simpson, Andrew Gelman
August 2018
1 Introduction
\(\widehat{R}\) (Gelman and Rubin, 1992; Brooks and Gelman, 1998) and split-\(\widehat{R}\) (Gelman et al., 2013) are based on the ratio of between and within-sequence marginal variances of the simulations. Effective sample size estimate (\(N_{\mbox{eff}}\)) by Gelman et al. (2013) is based on marginal variance and autocorrelation estimates from split chains. If a marginal distribution of some parameter does not have finite mean and variance, \(\widehat{R}\), split-\(\widehat{R}\) and \(N_{\mbox{eff}}\) for that parameter are not well defined.
To have a well defined MCMC convergence diagnostic and generic relative MCMC efficiency estimate, we propose to compute split-\(\widehat{R}\) and \(R_{\mbox{eff}}=N_{\mbox{eff}}/N\) for rank normalized simulations. For marginals with finite mean and variance, Monte Carlo error estimates for the desired expectations can be computed as before. In addition we propose two local efficiency measures based on relative efficiency of quantiles and small interval probabilities.
1.1 Rank normalization
Rank normalization is commonly used in Spearman’s rank correlation coefficient, but also quite commonly used to transform predictor values in epidemiology.
Rank normalization computes ranks for draws from all chains combined, applies z-transformation to ranks, and reconstructs the chains. Rank normalization makes the subsequent computations well defined and bijective invariant. This means that we get same results, for example, if we use unconstrained or constrained parameter. Furthermore to improve sensitivity to scale differences between chains, we consider folded rank-normalization which uses rank distance from the median.
1.2 Relative efficiency of median and other quantile estimates
Classic \(N_{\mbox{eff}}\) focused on posterior mean, but as it can be non-finite, we may prefer to report median and other quantiles by default. If median is \(x'\), then \(p(x < x')=0.5\). We can’t present this as a simple expectation, but probability \(p(x < x')\) can be estimated as \(\frac{1}{S}\sum I_{x<x'}\). The indicator function transforms simulation draws to 0’s and 1’s, and thus the subsequent computations are bijective invariant. Instead of estimating the efficiency of median directly, we can estimate the efficiency of \(p(x < \hat{Q}_{0.5})\), where \(\hat{Q}_{0.5}\) is the empirical median from combined draws. This should be good approximation, if the efficiency for \(p(x < \hat{Q}_{\alpha})\) is smoothly changing when \(\alpha\) is close to \(0.5\).
We can compute approximative relative efficiency for median absolute deviation based scale estimate (MAD_SD), by computing the relative efficiency for median of absolute deviations from median of all draws.
We can compute relative efficiency for other quantiles than median. We observe that often the relative efficiency for tail quantiles is lower than for mean or median.
1.3 Relative efficiency of small interval probability estimates
Next examine the relative efficiency in different parts of the posterior by computing the relative efficiency for \(p(\hat{Q}_\alpha \le x < \hat{Q}_{\alpha+\delta})\), where \(\hat{Q}_\alpha\) is empirical \(\alpha\)-quantile, \(delta=1/M\) is the length of the interval with some positive integer \(M\), and \(\alpha \in (0,\delta,\ldots,1-\delta)\). Each interval has \(N/M\) draws, and the efficiency measures autocorrelation of an indicator function which indicates when the values are inside the specific interval. This gives us local efficiency measure which does not depend on the shape of the distribution. We demonstrate it’s usefulness as a diagnostic tool and educational tool to teach that why expectations of different functionals may have different \(N_{\mbox{eff}}\).
1.4 Rank plots
We illustrate that rank plots can be used as more compact alternative to trace plots.
1.5 Finite mean and variance of a marginal
There will be a separate case study for diagnostics to check when certain marginal has finite mean and variance needed in order to report posterior mean and MCSE.
1.6 Abbreviations
In the code and figures the following abbreviations have been used
- Rhat = classic no-split-Rhat
- sRhat = classic split-Rhat
- zsRhat = rank-normalized split-Rhat
- all draws are jointly ranked, z-transformed and put back to the chains
- can detect differences in location and trends
- fzsRhat = folded rank-normalized split-Rhat
- all draws are jointly “folded” by computing absolute deviation from median, ranked, z-transformed and put back to the chains
- can detect differences in scales
- neff = no-split effective sample size
- reff = neff / N
- zsneff = rank-normalized split effective sample size
- estimates the efficiency of mean estimate for rank normalized values
- zsreff = zsneff / N
- fzsneff = folded rank-normalized split effective sample size
- estimates the efficiency of rank normalized mean absolute deviation
- fzsreff = fzsneff / N
- medsneff = median split effective sample size
- estimates the efficiency of median
- approximated using efficiency for \(p(x < \hat{Q}_{0.5})\)
- medsreff = medsneff / N
- madsneff = mad split effective sample size
- estimates the efficiency of median absolute deviation
- madsreff = madsneff / N
1.7 Proposed change in Stan
The proposal is to switch in Stan
- from split-\(\widehat{R}\) (sRhat) to maximum of rank-normalized-split-\(\widehat{R}\) and folded-rank-normalized-split-\(\widehat{R}\) (max(zsRhat, fzsRhat))
- from classic effective sample size estimate (neff, which currently doesn’t use chain splitting) to rank-normalized-split-\(R_{\mbox{eff}}\) (zsreff) and folded-rank-normalized-split-\(R_{\mbox{eff}}\) (fzsreff)
- if computing median and MAD_SD, report corresponding \(R_{\mbox{eff}}\)’s (medsreff, madsreff)
Justifications for the change are
Rank normalization makes zsRhat and zsreff well defined for all distributions, bijective transformation invariant and more stable than sRhat and neff. Adding folded version fzsRhat helps to detect differences in scale.
CDF based medsreff and masreff are well defined for all distributions, bijective transformation invariant, more stable than neff, and focuses on the median and mad_sd we are reporting.
In displays we propose eventually to use Rhat to denote the new version, but to make it more clear that rank-normalized efficiency is different we would switch to use R_eff and display relative efficiency. The relative efficiency is also easier to check for the low values as we don’t need to check the total sample size.
1.8 Warning thresholds
Based on the experiments, more strict convergence diagnostic and relative efficiency warning limits could be used.
We propose following warning thresholds, although additional experiments would be useful:
- 1.01 or 1.02 for new Rhat max(zsRhat, fzsRhat)
- 0.1 for new relative efficiency (zsreff, fzsreff, medsreff, madsreff)
Plots have dashed lines based on these thresholds (in continuous plots, a dashed line at 1.005 is plotted instead of 1.01, as values larger than that are usually rounded in our summaries to 1.01).
2 Setup
Load libraries
library(rstan)
rstan_options(auto_write = TRUE)
library(bayesplot)
library(grid)
library(gridExtra)
source('monitornew.R')
blank <- grid.rect(gp=gpar(col="white"))
3 Independent draws from normal distribution with additional trend, shift or scaling
This part focuses on diagnostics for
- all chains have a trend and similar marginal distribution
- one of the chains has different mean
- one of the chains has lower marginal variance
To simplify, in this part independent draws are used as a proxy for very efficient MCMC. In addition the draws are from normal distribution and we discuss the behavior for non-Gaussian distributions later.
3.1 Trend: All chains from the same \(N(0,1)\) distribution plus linear trend added to all chains.
res <- data.frame(iters=integer(),trend=double(),rep=integer(),sRhat=double(),zsRhat=double(),fzsRhat=double(),neff=integer(),zsneff=integer())
chains=4
for (iters in c(250, 1000, 4000)) {
for (a in c(0, 0.25, 0.5, 0.75, 1)) {
for (irep in 1:10) {
r <- array(rnorm(iters*chains), c(iters,chains))
r <- r + seq(-a, a, length.out=iters)
rs <- monitornew(r, warmup=0, digits_summary = 2, probs=c(0.5), print = FALSE)
res<-rbind(res,data.frame(iters=iters,trend=a,rep=irep,Rhat=rs[,'Rhat'],sRhat=rs[,'sRhat'],zsRhat=rs[,'zsRhat'],fzsRhat=rs[,'fzsRhat'],neff=rs[,'neff'],zsneff=rs[,'zsneff'],reff=rs[,'reff'],zsreff=rs[,'zsreff']))
}
}
}If we don’t split chains, Rhat misses trends if all chains still have similar marginal distribution.
ggplot(data=res, aes(y=Rhat, x=trend)) + geom_point() + geom_jitter() + facet_grid(. ~ iters) + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1) + ggtitle('Rhat without splitting chains')
Split-Rhat can detect trends, even if the marginals of chains are similar.
ggplot(data=res, aes(y=zsRhat, x=trend)) + geom_point() + geom_jitter() + facet_grid(. ~ iters) + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1) + ggtitle('Split-Rhat')
Result: Split-Rhat is useful for detecting non-stationarity in chains. If we use threshold 1.01, we can detect trends which have variance 2% from the marginal variance. If we use threshold 1.1, we can detect trends which have variance 30% from the marginal variance.
Relative efficiency (effective sample size divided by the number of draws) is based on split Rhat and autocorrelations.
ggplot(data=res, aes(y=zsreff, x=trend)) + geom_point() + geom_jitter() + facet_grid(. ~ iters) + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + ggtitle('Relative efficiency (zsreff)') + scale_y_continuous(breaks = seq(0,1,by=0.25))
Result: Split-Rhat is more sensitive to trends in small sample sizes, but relative efficiency is more sensitive with larger samples sizes (as autocorrelations can be estimated more accurately).
Advice: If in doubt, run longer chains for more accurate estimates.
3.2 Shifted: Three chains from \(N(0,1)\) and one chain from \(N(a,1)\)
Next we investigate the sensitivity to detect if one of the chains has not converged to the same distribution as the others, but has a different mean.
res <- data.frame(iters=integer(),trend=double(),rep=integer(),sRhat=double(),zsRhat=double(),fzsRhat=double(),neff=integer(),zsneff=integer())
chains=4
for (iters in c(250, 1000, 4000)) {
for (a in c(0, 0.25, 0.5, 0.75, 1)) {
for (irep in 1:10) {
r <- array(rnorm(iters*chains), c(iters,chains))
r[,1] <- r[,1] + a
rs <- monitornew(r, warmup=0, digits_summary = 2, probs=c(0.5), print = FALSE)
res<-rbind(res,data.frame(iters=iters,shift=a,rep=irep,Rhat=rs[,'Rhat'],sRhat=rs[,'sRhat'],zsRhat=rs[,'zsRhat'],fzsRhat=rs[,'fzsRhat'],neff=rs[,'neff'],zsneff=rs[,'zsneff'],reff=rs[,'reff'],zsreff=rs[,'zsreff']))
}
}
}ggplot(data=res, aes(y=zsRhat, x=shift)) + geom_point() + geom_jitter() + facet_grid(. ~ iters) + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1) + ggtitle('Split-Rhat')
Result: If we use threshold 1.01, we can detect shift with magnitude one third of marginal standard deviation. If we use threshold 1.1, we can detect shift with the magnitude equal to marginal standard deviation. The rank plot can be used to visualize where the problem is.
ggplot(data=res, aes(y=zsreff, x=shift)) + geom_point() + geom_jitter() + facet_grid(. ~ iters) + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + ggtitle('Relative efficiency (zsreff)') + scale_y_continuous(breaks = seq(0,1,by=0.25))
Result: The relative efficiency is not as sensitive, but a shift with magnitude half of the marginal standard deviation will lead to very low efficiency when sample size increases.
Rank plot visualisation for a case with 4 chains and 250 draws per chain, and shift = 0.5.
iters=250;chains=4;a=0.5;
r <- array(rnorm(iters*chains), c(iters,chains))
r[,1] <- r[,1] + a
colnames(r)<-c('1','2','3','4');
mcmc_hist(r_scale(r), breaks = seq(0,iters*chains,by=iters*chains/20)+0.5)
Rhat is less than 1.05, but the rank plot clearly shows where the problem is.
3.3 Scaled: Three chains from \(N(0,1)\) and one chain from \(N(0,a)\)*
Next we investigate the sensitivity to detect if one of the chains has not converged to the same distribution as the others but has lower marginal variance.
res <- data.frame(iters=integer(),trend=double(),rep=integer(),sRhat=double(),zsRhat=double(),fzsRhat=double(),neff=integer(),zsneff=integer())
chains=4
for (iters in c(250, 1000, 4000)) {
for (a in c(0, 0.25, 0.5, 0.75, 1)) {
for (irep in 1:10) {
r <- array(rnorm(iters*chains), c(iters,chains))
r[,1] <- r[,1]*a
rs <- monitornew(r, warmup=0, digits_summary = 2, probs=c(0.5), print = FALSE)
res<-rbind(res,data.frame(iters=iters,scale=a,rep=irep,Rhat=rs[,'Rhat'],sRhat=rs[,'sRhat'],zsRhat=rs[,'zsRhat'],fzsRhat=rs[,'fzsRhat'],neff=rs[,'neff'],zsneff=rs[,'zsneff'],reff=rs[,'reff'],zsreff=rs[,'zsreff'],fzsreff=rs[,'fzsreff']))
}
}
}ggplot(data=res, aes(y=zsRhat, x=scale)) + geom_point() + geom_jitter() + facet_grid(. ~ iters) + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1) + ggtitle('Split-Rhat')
Result: Split-Rhat is not able to detect differences in scale.
ggplot(data=res, aes(y=fzsRhat, x=scale)) + geom_point() + geom_jitter() + facet_grid(. ~ iters) + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1) + ggtitle('Folded-split-Rhat')
Result: Folded-Split-Rhat focuses instead to scales, and can detect scale differences.
Result: If we use threshold 1.01, we can detect a chain with scale less than 3/4 of the scale of the others. If we use threshold 1.1, we can detect a chain with scale less than 1/4 of the scale of the others.
ggplot(data=res, aes(y=zsreff, x=scale)) + geom_point() + geom_jitter() + facet_grid(. ~ iters) + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + ggtitle('Relative efficiency (zsreff)') + scale_y_continuous(breaks = seq(0,1,by=0.25))
Result: The relative efficiency does not see a problem as it focuses on location.
ggplot(data=res, aes(y=fzsreff, x=scale)) + geom_point() + geom_jitter() + facet_grid(. ~ iters) + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + ggtitle('Folded relative efficiency (fzsreff)') + scale_y_continuous(breaks = seq(0,1,by=0.25))
Result: The relative efficiency of sd does see the problem as it focuses on scale .
Rank plot visualisation for a case with 4 chains and 250 draws per chain, and one chains has scale 0.75.
iters=250;chains=4;a=0.75;
r <- array(rnorm(iters*chains), c(iters,chains))
r[,1] <- r[,1]*a
colnames(r)<-c('1','2','3','4');
mcmc_hist(r_scale(r), breaks = seq(0,iters*chains,by=iters*chains/20)+0.5)
Folded Rhat is 1.06, but the rank plot clearly shows where the problem is.
4 Cauchy: mean and variance not finite
Classic Rhat is based on calculating within and between chain variances. If the marginal distribution of the chain is such that variance is not defined classic Rhat is not well justified. Next we use Cauchy as an example of such distribution. Also in cases where mean and variance are finite, the distribution can be far from Gaussian and especially distributions with very long tails cause instability for variance and autocorrelation estimates. To alleviate these problems we Split-Rhat for rank-normalized draws. Ranks are computed by using average for ties, as it is deterministic and constant chains are mapped to constant chains.
The following Cauchy models are from Michael Betancourt’s case study Fitting The Cauchy Distribution
4.1 Nominal parameterization of Cauchy
The the nominal model with direct parameterization is as follows.
writeLines(readLines("cauchy_nom.stan"))parameters {
vector[50] x;
}
model {
x ~ cauchy(0, 1);
}
generated quantities {
real I = fabs(x[1]) < 1 ? 1 : 0;
}4.1.1 Default Stan options
Run the nominal model
fit_nom <- stan(file='cauchy_nom.stan', seed=7878, refresh = 0)Warning: There were 1749 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 10. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceededWarning: There were 2 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-lowWarning: Examine the pairs() plot to diagnose sampling problemsThe diagnostics warn about very large number of transitions after warmup that exceeded the maximum treedepth, which is likely due to very long tail of Cauchy. All chains have low estimated Bayesian Fraction of Missing Information indicating also slow mixing.
MCMC trace for the first parameter looks wild with occasional large values
samp<-as.array(fit_nom)
mcmc_trace(samp[,,1])
Let’s check Rhat and relative efficiency diagnostics.
res<-monitornew(samp[,,1:50], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
p1 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`sRhat`)) + geom_point() + ggtitle('Classic Split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1) + ylim(c(.99 ,1.26))
p2 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`zsRhat`)) + geom_point() + ggtitle('Rank normalized split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1) + ylim(c(.99,1.26))
p3 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`fzsRhat`)) + geom_point() + ggtitle('Folded rank split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1) + ylim(c(.99,1.26))
grid.arrange(p1, p2, p3, nrow = 1)
For one parameter Rhats exceed the classic threshold 1.1. Depending on Rhat, a few others exceed the threshold 1.01. The interpretation of classic Rhat is however not well defined. The rank normalized split-Rhat have several values over 1.01.
res<-monitornew(samp[,,1:50], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
p1 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=reff)) + geom_point() + ggtitle('Classic rel. efficiency (reff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(breaks = seq(0,1,by=0.25), limits = c(0,1.1))
p2 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=zsreff)) + geom_point() + ggtitle('New rel. efficiency (zsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(breaks = seq(0,1,by=0.25), limits = c(0,1.1))
p3 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=fzsreff)) + geom_point() + ggtitle('Folded rel. efficiency (fzsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(breaks = seq(0,1,by=0.25), limits = c(0,1.1))
p4 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=medsreff)) + geom_point() + ggtitle('Rel.eff. of median (medsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(breaks = seq(0,1,by=0.25), limits = c(0,1.1))
p5 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=madsreff)) + geom_point() + ggtitle('Rel.eff. of mad (madsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(breaks = seq(0,1,by=0.25), limits = c(0,1.1))
grid.arrange(p1, p2, p3, blank, p4, p5, nrow = 2)
Both classic and rank normalized relative efficiency estimates have several near zero values, and the overall samples shouldn’t be trusted.
Result: Relative efficiency is more sensitive than (rank-normalized) split-Rhat to detect problems of slow mixing.
We check the indicator function = and lp__
res<-monitornew(samp[,,51:52], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
cat('I: r_eff = ', round(res['I','zsreff'],2), '\n')I: r_eff = 0.14 cat('lp__: r_eff = ', round(res['lp__','zsreff'],2), '\n')lp__: r_eff = 0.05 Relative efficiency for lp__ is worryingly low.
We examine the relative efficiency in different parts of the posterior, by computing the relative efficiency for \(p(\hat{Q}_\alpha \le x < \hat{Q}_{\alpha+0.05})\), where \(\hat{Q}_\alpha\) is empirical \(\alpha\)-quantile and \(\alpha \in (0,0.05,\ldots,0.95)\). Each interval has \(5%\) of draws, and the efficiency measures autocorrelation of an indicator function which indicates when the values are inside the specific interval. This gives us local efficiency measure which is does not depend on the shape of the distribution.
plot_local_reff <- function(samp = NULL, parami = NULL, nalpha=NULL) {
zsreffs <- double()
delta <- 1/nalpha
for (alpha in seq(0, 1-delta, by=delta)) {
rs <- monitornew(samp[,,parami]>=quantile(samp[,,parami],alpha)&samp[,,parami]<quantile(samp[,,parami],alpha+delta),warmup = 0, probs = 0.5, digits_summary = 2, print = FALSE)
zsreffs <- c(zsreffs, rs[1,'zsreff'])
}
df <- data.frame(quantile=seq(0, 1, by=delta), zsreff=c(zsreffs,zsreffs[nalpha]))
ggplot(data=df, aes(x=quantile, y=zsreff)) +
geom_step() +
geom_hline(yintercept = c(0,1)) +
geom_hline(yintercept = 0.1, linetype = 'dashed') +
scale_x_continuous(breaks=seq(0,1, by=0.1)) +
scale_y_continuous(breaks = seq(0,1,by=0.25), limits = c(0,1.1))
}
plot_local_reff(samp = samp, parami = 1, nalpha = 20)
We see that MCMC is less efficient in tails.
Alternative way to examine the relative efficiency in different parts of the posterior, is to compute relative efficiency for \(p(x < \hat{Q}_{\alpha})\), where \(\hat{Q}_\alpha\) is empirical \(\alpha\)-quantile and \(\alpha \in (0.025,\ldots,0.975)\). This is related to the efficiency of quantile estimates. Each interval has specified proportion of draws, and the efficiency measures autocorrelation of an indicator function which indicates when the values are inside the specific interval. This gives us local efficiency measure which is does not depend on the shape of the distribution.
plot_quantile_reff <- function(samp = NULL, parami = NULL, nalpha=NULL) {
zsreffs <- double()
delta <- 1/nalpha
for (alpha in seq(delta, 1-delta, by=delta)) {
rs <- monitornew(samp[,,parami]<quantile(samp[,,parami],alpha),warmup = 0, probs = 0.5, digits_summary = 2, print = FALSE)
zsreffs <- c(zsreffs, rs[1,'zsreff'])
}
df <- data.frame(quantile=seq(delta, 1-delta, by=delta), zsreff=zsreffs)
ymax<-max(1,round(max(zsreffs)+0.15,1))
ggplot(data=df, aes(x=quantile, y=zsreff)) +
geom_point() +
geom_hline(yintercept = c(0,1)) +
geom_hline(yintercept = 0.1, linetype = 'dashed') +
scale_x_continuous(breaks=seq(0,1, by=0.1)) +
scale_y_continuous(breaks = seq(0,ymax,by=0.25), limits = c(0,ymax))
}
plot_quantile_reff(samp = samp, parami = 1, nalpha = 40)
Rank plot visualisation of x[1], which has relative efficiency 0.39
mcmc_hist(r_scale(samp[,,1]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/50)+0.5)
Chains have clearly different rank plots.
Rank plot visualisation of x[7], which has relative efficiency 0.02
mcmc_hist(r_scale(samp[,,7]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/50)+0.5)
There are clear problems with mixing.
4.1.2 Default Stan options + max_treedepth=20
We can try to improve the performance by increasing the maximum treedepth.
fit_nom_td20 <- stan(file='cauchy_nom.stan', seed=7878, refresh = 0, control = list(max_treedepth=20))Warning: There were 3 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-lowWarning: Examine the pairs() plot to diagnose sampling problemsMCMC trace for the first parameter looks wild with occasional large values
samp<-as.array(fit_nom_td20)
mcmc_trace(samp[,,1])
Treedepth exceedence and Bayesian Fraction of Missing Information are HMC/NUTS specific diagnostics. Let’s see what generic Split-Rhat and relative efficiency can detect.
Check the diagnostics for all \(x\)’s.
res<-monitornew(samp[,,1:50], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
p1 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`sRhat`)) + geom_point() + ggtitle('Classic Split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1) + ylim(0.999, 1.045)
p2 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`zsRhat`)) + geom_point() + ggtitle('Rank normalized split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1) + ylim(0.999, 1.045)
p3 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`fzsRhat`)) + geom_point() + ggtitle('Folded rank split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1) + ylim(0.999, 1.045)
grid.arrange(p1, p2, p3, nrow = 1)
All Rhats are below 1.1, but many are over 1.01. Classic Rhat has more variation than the rank normalized (and the classic is not well defined). The folded rank normalized Rhat shows that there is still more variation in scale than location between different chains.
res<-monitornew(samp[,,1:50], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
p1 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=reff)) + geom_point() + ggtitle('Classic rel. efficiency (reff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(breaks = seq(0,1.5,by=0.25), limits = c(0,1.1))
p2 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=zsreff)) + geom_point() + ggtitle('New rel. efficiency (zsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(breaks = seq(0,1.5,by=0.25), limits = c(0,1.1))
p3 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=fzsreff)) + geom_point() + ggtitle('New rel. efficiency (zsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(breaks = seq(0,1.5,by=0.25), limits = c(0,1.1))
p4 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=medsreff)) + geom_point() + ggtitle('Rel.eff. of median (medsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(breaks = seq(0,1,by=0.25), limits = c(0,1.1))
p5 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=madsreff)) + geom_point() + ggtitle('Rel.eff. of mad (madsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(breaks = seq(0,1,by=0.25), limits = c(0,1.1))
grid.arrange(p1, p2, p3, blank, p4, p5, nrow = 2)Warning: Removed 1 rows containing missing values (geom_point).
Some classic relative efficiencies (reffs) are close to zero. If we wouldn’t realize that the variance is infinite, we might try to run longer chains, but in case of infinite variance, zero efficiency is the truth and longer chains don’t help. However other quantities can be well defined, and that’s why it is useful also to look at the rank normalized version as a generic transformation with finite mean and variance. The smallest rank-normalized r_eff are around 0.25, which is not that bad. There is still large variation between parameters although they all have the same distribution.
Result: Rank normalized relative efficiency is more stable than not-rank-normalized which is not well defined for Cauchy.
We check the indicator function = and lp__
res<-monitornew(samp[,,51:52], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
cat('I: r_eff = ', round(res['I','zsreff'],2), '\n')I: r_eff = 0.13 cat('lp__: r_eff = ', round(res['lp__','zsreff'],2), '\n')lp__: r_eff = 0.06 Although increasing max treedepth improved diagnostics for x, the efficiency for I and lp__ didn’t change.
We examine the relative efficiency in different parts of the posterior, by computing the relative efficiency for \(p(\hat{Q}_\alpha \le x < \hat{Q}_{\alpha+0.05})\), where \(\hat{Q}_\alpha\) is \(\alpha\)-quantile and \(\alpha \in (0,0.05,\ldots,0.95)\).
plot_local_reff(samp = samp, parami = 1, nalpha = 20)
Increasing max treedepth improved the efficiency slightly in the tails.
We examine the relative efficiency in different parts of the posterior, by computing the relative efficiency for \(p(x < \hat{Q}_{\alpha+0.05})\), where \(\hat{Q}_\alpha\) is \(\alpha\)-quantile and \(\alpha \in (0,0.025,\ldots,0.975)\).
plot_quantile_reff(samp = samp, parami = 1, nalpha = 40)
Rank plot visualisation of x[38], which has the smallest relative efficiency 0.23 among x’s.
mcmc_hist(r_scale(samp[,,38]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/50)+0.5)
Clear problems.
Rank plot visualisation of lp__, which has relative efficiency 0.06.
mcmc_hist(r_scale(samp[,,52]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/50)+0.5)
Doesn’t look so good.
4.1.3 Default Stan options + max_treedepth=20 + 8 times longer sampling
Let’s try running longer chain
fit_nom_td20l <- stan(file='cauchy_nom.stan', seed=7878, refresh = 0, control = list(max_treedepth=20), warmup=1000, iter = 9000)Warning: There were 3 transitions after warmup that exceeded the maximum treedepth. Increase max_treedepth above 20. See
http://mc-stan.org/misc/warnings.html#maximum-treedepth-exceededWarning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-lowWarning: Examine the pairs() plot to diagnose sampling problemsMCMC trace for the first parameter looks wild with occasional large values
samp<-as.array(fit_nom_td20l)
mcmc_trace(samp[,,1])
Treedepth exceedence and Bayesian Fraction of Missing Information are HMC/NUTS specific diagnostics. Let’s see what generic Split-Rhat and relative efficiency can detect. Check the diagnostics for all \(x\)’s.
res<-monitornew(samp[,,1:50], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
p1 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`sRhat`)) + geom_point() + ggtitle('Classic split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1) + ylim(.999, 1.009)
p2 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`zsRhat`)) + geom_point() + ggtitle('Rank norm. split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1) + ylim(.999, 1.009)
p3 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`fzsRhat`)) + geom_point() + ggtitle('Fold. rank split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1) + ylim(.999, 1.009)
grid.arrange(p1, p2, p3, nrow = 1)
All Rhats are below 1.01. Classic Rhat has more variation than the rank normalized (but the classic is not well defined).
res<-monitornew(samp[,,1:50], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
p1 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=reff)) + geom_point() + ggtitle('Classic rel. efficiency (reff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + ylim(0,1.1)
p2 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=zsreff)) + geom_point() + ggtitle('New rel. efficiency (zsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + ylim(0,1.1)
p3 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=fzsreff)) + geom_point() + ggtitle('New rel. efficiency (zsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + ylim(0,1.1)
p4 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=medsreff)) + geom_point() + ggtitle('Rel.eff. of median (medsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(breaks = seq(0,1,by=0.25), limits = c(0,1.1))
p5 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=madsreff)) + geom_point() + ggtitle('Rel.eff. of mad (madsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(breaks = seq(0,1,by=0.25), limits = c(0,1.1))
grid.arrange(p1, p2, p3, blank, p4, p5, nrow = 2)
Most classic relative efficiencies (reffs) are close to zero. If we wouldn’t realize that the variance is infinite, we might try to run longer chains, but in case of infinite variance zero is the true value and longer chains don’t help. However other quantities can be well defined, and that’s why it can be useful also to look at the rank normalized version as a generic transformation with finite mean and variance. The rank-normalized r_eff’s start to concentrate around 0.6, which is quite good. There is still some variation between parameters.
Result: Rank normalized relative efficiency is more stable than not-rank-normalized which is not well defined for Cauchy.
We check the indicator function = and lp__
res<-monitornew(samp[,,51:52], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
cat('I: r_eff = ', round(res['I','reff'],2), '\n')I: r_eff = 0.16 cat('lp__: r_eff = ', round(res['lp__','zsreff'],2), '\n')lp__: r_eff = 0.04 We examine the relative efficiency in different parts of the posterior, by computing the relative efficiency for \(p(\hat{Q}_\alpha \le x < \hat{Q}_{\alpha+0.05})\), where \(\hat{Q}_\alpha\) is \(\alpha\)-quantile and \(\alpha \in (0,0.05,\ldots,0.95)\).
plot_local_reff(samp = samp, parami = 1, nalpha = 20)
Increasing chain length did not seem to change the relative efficiency. With more draws from the longer chains we can use finer resolution for the local efficiency estimates.
plot_local_reff(samp = samp, parami = 1, nalpha = 100)
The efficiency far in the tails is worryingly low.
The alternative local relative efficiency.
plot_quantile_reff(samp = samp, parami = 1, nalpha = 100)
Rank plot visualisation of x[34], which has the smallest relative efficiency 0.22 among x’s.
mcmc_hist(r_scale(samp[,,34]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/100)+0.5)
Still problems at the edges.
Rank plot visualisation of lp__, which has relative efficiency 0.04.
mcmc_hist(r_scale(samp[,,52]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/100)+0.5)
There seems to be slight problems also in mixing for energy.
4.2 1st alternative parameterization of Cauchy
Next we examine alternative parameterization and consider Cauchy as a scale mixture of Gaussian distributions. The model has two parameters and the Cauchy distributed x’s can be computed from those. In addition of improved sampling performance, the example illustrates that focus of diagnostics matter.
writeLines(readLines("cauchy_alt_1.stan"))parameters {
vector[50] x_a;
vector<lower=0>[50] x_b;
}
transformed parameters {
vector[50] x = x_a ./ sqrt(x_b);
}
model {
x_a ~ normal(0, 1);
x_b ~ gamma(0.5, 0.5);
}
generated quantities {
real I = fabs(x[1]) < 1 ? 1 : 0;
}Run the alternative model
fit_1 <- stan(file='cauchy_alt_1.stan', seed=7878, refresh = 0)There are no warnings, and the sampling is much faster.
samp<-as.array(fit_1)
res<-monitornew(samp[,,101:150], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
p1 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`sRhat`)) + geom_point() + ggtitle('Classic split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1)
p2 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`zsRhat`)) + geom_point() + ggtitle('Rank norm. split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1)
p3 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`fzsRhat`)) + geom_point() + ggtitle('Fold. rank split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1)
grid.arrange(p1, p2, p3, nrow = 1)
All Rhats are below 1.01. Classic Rhat’s look also good even if the classic is not well defined for Cauchy distributed.
res<-monitornew(samp[,,101:150], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
p1 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=reff)) + geom_point() + ggtitle('Classic rel. efficiency (reff)') + geom_hline(yintercept = c(0,1)) + ylim(0,3.4)
p2 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=zsreff)) + geom_point() + ggtitle('New rel. efficiency (zsreff)') + geom_hline(yintercept = c(0,1)) + ylim(0,3.4)
p3 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=fzsreff)) + geom_point() + ggtitle('New rel. efficiency (zsreff)') + geom_hline(yintercept = c(0,1)) + ylim(0,3.4)
p4 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=medsreff)) + geom_point() + ggtitle('Rel.eff. of median (medsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(limits = c(0,3.4))
p5 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=madsreff)) + geom_point() + ggtitle('Rel.eff. of mad (madsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(limits = c(0,3.4))
grid.arrange(p1, p2, p3, blank, p4, p5, nrow = 2)
Result: Rank normalized r_eff’s have less variation than classic one which is not well defined for Cauchy.
We check the indicator function = and lp__
res<-monitornew(samp[,,151:152], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
cat('I: r_eff = ', round(res['I','zsreff'],2), '\n')I: r_eff = 0.75 cat('lp__: r_eff = ', round(res['lp__','zsreff'],2), '\n')lp__: r_eff = 0.33 The relative efficiency for these are also much better than with nominal parameterization.
We examine the relative efficiency in different parts of the posterior, by computing the relative efficiency for \(p(\hat{Q}_\alpha \le x < \hat{Q}_{\alpha+0.05})\), where \(\hat{Q}_\alpha\) is \(\alpha\)-quantile and \(\alpha \in (0,0.05,\ldots,0.95)\).
samp <- as.array(fit_1)
plot_local_reff(samp = samp, parami = 101, nalpha = 20)
The relative efficiency is good in all parts of the posterior.
We examine the relative efficiency in different parts of the posterior, by computing the relative efficiency for \(p(x < \hat{Q}_{\alpha+0.05})\), where \(\hat{Q}_\alpha\) is \(\alpha\)-quantile and \(\alpha \in (0,0.025,\ldots,0.975)\).
samp <- as.array(fit_1)
plot_quantile_reff(samp = samp, parami = 101, nalpha = 40)
res<-monitornew(samp[,,101:150], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
res1<-monitornew(samp[,,1:50], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
res2<-monitornew(samp[,,51:100], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
cat('Mean r_eff for xs = ' , round(mean(res[,'zsreff']),2), '\n')Mean r_eff for xs = 1.03 cat('Mean fr_eff for xs = ' , round(mean(res[,'fzsreff']),2), '\n')Mean fr_eff for xs = 0.52 cat('Mean r_eff for x_as = ' , round(mean(res1[,'zsreff']),2), '\n')Mean r_eff for x_as = 1.15 cat('Mean r_eff for x_bs = ' , round(mean(res2[,'zsreff']),2), '\n')Mean r_eff for x_bs = 0.69 Result: We see that, the relative efficiency of the interesting x, can be different from the relative efficiency of the parameters that used to compute it. Also r_eff we are computing is for rank normalized and to compute Monte Carlo errors, the final target function should be taken into account, but we just want to avoid the complication of the potentially non-existing moments at this point, and come back later in other paper how to use Pareto diagnostic to check whether necessary moments exist to compute the desired expectations and their Monte Carlo errors.
Rank plot visualisation of x[9], which has the smallest relative efficiency 0.85 among x’s.
mcmc_hist(r_scale(samp[,,109]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/50)+0.5)
Looks better than for the nominal parameterization.
Rank plot visualisation of lp__, which has relative efficiency 0.33.
mcmc_hist(r_scale(samp[,,152]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/50)+0.5)
Looks better than for the nominal parameterization. It would be useful to have the reference plotted.
4.3 Another alternative parameterization of Cauchy
Another alternative parameterization is univariate transformation as shown in the following code.
writeLines(readLines("cauchy_alt_3.stan"))parameters {
vector<lower=0, upper=1>[50] x_tilde;
}
transformed parameters {
vector[50] x = tan(pi() * (x_tilde - 0.5));
}
model {
// Implicit uniform prior on x_tilde
}
generated quantities {
real I = fabs(x[1]) < 1 ? 1 : 0;
}4.4 2nd alternative parameterization of Cauchy
Run the alternative model (3rd alternative in Michael Betancourt’s case study)
fit_3 <- stan(file='cauchy_alt_3.stan', seed=7878, refresh = 0)There are no warnings, and the sampling is much faster than for the nominal model.
samp<-as.array(fit_3)
res<-monitornew(samp[,,51:100], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
p1 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`sRhat`)) + geom_point() + ggtitle('Classic split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1) + ylim(0.999, 1.007)
p2 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`zsRhat`)) + geom_point() + ggtitle('Rank norm. split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1) + ylim(0.999, 1.007)
p3 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`fzsRhat`)) + geom_point() + ggtitle('Fold. rank split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1) + ylim(0.999, 1.007)
grid.arrange(p1, p2, p3, nrow = 1)
All Rhats except some folded Rhats are below 1.01. Classic Rhat’s look also good even if the classic is not well defined for Cauchy distributed.
res<-monitornew(samp[,,51:100], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
p1 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`reff`)) + geom_point() + ggtitle('Classic rel. efficiency (reff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + ylim(0,3.4)
p2 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`zsreff`)) + geom_point() + ggtitle('New rel. efficiency (zsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + ylim(0,3.4)
p3 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`fzsreff`)) + geom_point() + ggtitle('New rel. efficiency (zsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + ylim(0,3.4)
p4 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=medsreff)) + geom_point() + ggtitle('Rel.eff. of median (medsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(limits = c(0,3.4))
p5 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=madsreff)) + geom_point() + ggtitle('Rel.eff. of mad (madsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(limits = c(0,3.4))
grid.arrange(p1, p2, p3, blank, p4, p5, nrow = 2)
Result: Rank normalized r_eff’s have less variation than classic one which is not well defined for Cauchy. Rank normalized r_eff’s are slightly larger than 1, which is possible for antithetic Markov chains which have negative correlation for odd lags.
We check the indicator function = and lp__
res<-monitornew(samp[,,101:102], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
cat('I: r_eff = ', round(res['I','zsreff'],2), '\n')I: r_eff = 0.43 cat('lp__: r_eff = ', round(res['lp__','zsreff'],2), '\n')lp__: r_eff = 0.36 The relative efficiency for these are also much better than with nominal parameterization.
We examine the relative efficiency in different parts of the posterior, by computing the relative efficiency for \(p(\hat{Q}_\alpha \le x < \hat{Q}_{\alpha+0.05})\), where \(\hat{Q}_\alpha\) is \(\alpha\)-quantile and \(\alpha \in (0,0.05,\ldots,0.95)\).
samp <- as.array(fit_3)
plot_local_reff(samp = samp, parami = 51, nalpha = 20)
We examine the relative efficiency in different parts of the posterior, by computing the relative efficiency for \(p(x < \hat{Q}_{\alpha+0.05})\), where \(\hat{Q}_\alpha\) is \(\alpha\)-quantile and \(\alpha \in (0,0.025,\ldots,0.975)\).
samp <- as.array(fit_3)
plot_quantile_reff(samp = samp, parami = 51, nalpha = 40)
The relative efficiency in tails is worse than for the first alternative parameterization, although it’s still better than for the nominal model.
res<-monitornew(samp[,,51:100], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
res1<-monitornew(samp[,,1:50], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
cat('Mean r_eff for xs = ' , round(mean(res[,'zsreff']),2), '\n')Mean r_eff for xs = 1.25 cat('Mean fr_eff for xs = ' , round(mean(res[,'fzsreff']),2), '\n')Mean fr_eff for xs = 0.41 cat('Mean r_eff for x_as = ' , round(mean(res1[,'zsreff']),2), '\n')Mean r_eff for x_as = 1.25 Result: When the transformation is univariate and bijective the rank normalized Rhat and r_eff are transformation invariant.
Naturally when computing the desired expectations the final target function should be taken into account, but we just want to avoid the complication of the potentially non-existing moments at this point, and come back later in other paper how to use Pareto diagnostic to check whether necessary moments exist to compute the desired expectations and their Monte Carlo errors.
Rank plot visualisation of x[13], which has the smallest relative efficiency 1.0 among x’s.
mcmc_hist(r_scale(samp[,,63]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/50)+0.5)
Nothing special.
Rank plot visualisation of lp__, which has relative efficiency 0.33.
mcmc_hist(r_scale(samp[,,102]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/50)+0.5)
Nothing special.
4.5 Half-Cauchy with nominal parameterization
Half-Cauchy priors are common and, for example, in Stan usually set apparently using the nominal parameterization. However, when the constraint <lower=0> is used, Stan does the sampling automatically in the unconstrained log(x) space, which changes the geometry crucially.
writeLines(readLines("half_cauchy_nom.stan"))parameters {
vector<lower=0>[50] x;
}
model {
x ~ cauchy(0, 1);
}
generated quantities {
real I = fabs(x[1]) < 1 ? 1 : 0;
}Run the half-Cauchy with nominal parameterization (and positive constraint)
fit_half_nom <- stan(file='half_cauchy_nom.stan', seed=7878, refresh = 0)There are no warnings, and the sampling is much faster than for the Cauchy nominal model.
samp<-as.array(fit_half_nom)
res<-monitornew(samp[,,1:50], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
p1 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`sRhat`)) + geom_point() + ggtitle('Classic split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1)
p2 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`zsRhat`)) + geom_point() + ggtitle('Rank norm. split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1)
p3 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`fzsRhat`)) + geom_point() + ggtitle('Fold. rank split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1)
grid.arrange(p1, p2, p3, nrow = 1)
All Rhats are below 1.01. Classic Rhat’s look also good even if the classic is not well defined for half-Cauchy distribution.
res<-monitornew(samp[,,1:50], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
p1 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`reff`)) + geom_point() + ggtitle('Classic relative efficiency (reff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + ylim(0,2.2)
p2 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`zsreff`)) + geom_point() + ggtitle('New relative efficiency (zsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + ylim(0,2.2)
p3 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`fzsreff`)) + geom_point() + ggtitle('New relative efficiency (zsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + ylim(0,2.2)
p4 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=medsreff)) + geom_point() + ggtitle('Rel.eff. of median (medsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(limits = c(0,2.2))
p5 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=madsreff)) + geom_point() + ggtitle('Rel.eff. of mad (madsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(limits = c(0,2.2))
grid.arrange(p1, p2, p3, blank, p4, p5, nrow = 2)
Result: Rank normalized r_eff’s have less variation than classic one which is not well defined for Cauchy. Rank normalized r_eff’s are much larger than 1, which is possible for antithetic Markov chains which have negative correlation for odd lags.
Due to constraint ‘log(x) space, and we can also check the performance in that space.
res<-monitornew(log(samp[,,1:50]), warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
p1 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`reff`)) + geom_point() + ggtitle('Classic relative efficiency (reff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + ylim(0,2.2)
p2 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`zsreff`)) + geom_point() + ggtitle('New relative efficiency (zsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + ylim(0,2.2)
p3 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`fzsreff`)) + geom_point() + ggtitle('New relative efficiency (zsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + ylim(0,2.2)
p4 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=medsreff)) + geom_point() + ggtitle('Rel.eff. of median (medsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(limits = c(0,2.2))
p5 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=madsreff)) + geom_point() + ggtitle('Rel.eff. of mad (madsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(limits = c(0,2.2))
grid.arrange(p1, p2, p3, blank, p4, p5, nrow = 2)
\(log(x)\) is quite close to Gaussian, and thus classic r_eff is also close to rank normalized r_eff which is exactly same as for \(x\) as the rank normalized version is invariant to bijective transformations.
Result: Rank normalized r_eff is close to classic r_eff for transformations which make the distribution close to Gaussian.
We check the indicator function = and lp__
res<-monitornew(samp[,,51:52], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
cat('I: r_eff = ', round(res['I','zsreff'],2), '\n')I: r_eff = 1.33 cat('lp__: r_eff = ', round(res['lp__','zsreff'],2), '\n')lp__: r_eff = 0.27 Result: Relative efficiency for \(p(x_1<1)\) is larger than 1. This is possible for antithetic Markov chains, which have negative correlation for odd lags.
We examine the relative efficiency in different parts of the posterior, by computing the relative efficiency for \(p(\hat{Q}_\alpha \le x < \hat{Q}_{\alpha+0.05})\), where \(\hat{Q}_\alpha\) is \(\alpha\)-quantile and \(\alpha \in (0,0.05,\ldots,0.95)\).
samp <- as.array(fit_half_nom)
plot_local_reff(samp = samp, parami = 1, nalpha = 20)
The relative efficiency is good overall, with only a small dip in tails.
We examine the relative efficiency in different parts of the posterior, by computing the relative efficiency for \(p(x < \hat{Q}_{\alpha+0.05})\), where \(\hat{Q}_\alpha\) is \(\alpha\)-quantile and \(\alpha \in (0,0.025,\ldots,0.975)\).
samp <- as.array(fit_half_nom)
plot_quantile_reff(samp = samp, parami = 1, nalpha = 40)
Rank plot visualisation of x[32], which has the smallest relative efficiency 1.0 among x’s.
mcmc_hist(r_scale(samp[,,32]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/50)+0.5)
Looks like there are differences in scale of chains.
Rank plot visualisation of lp__, which has relative efficiency 0.33.
mcmc_hist(r_scale(samp[,,52]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/50)+0.5)
Maybe small differences in scales, but it’s difficult to know whether small variation from uniform is relevant.
4.6 Alternative parameterization of half-Cauchy
writeLines(readLines("half_cauchy_alt.stan"))parameters {
vector<lower=0>[50] x_a;
vector<lower=0>[50] x_b;
}
transformed parameters {
vector[50] x = x_a .* sqrt(x_b);
}
model {
x_a ~ normal(0, 1);
x_b ~ inv_gamma(0.5, 0.5);
}
generated quantities {
real I = fabs(x[1]) < 1 ? 1 : 0;
}Run half-Cauchy with alternative parameterization
fit_half_reparam <- stan(file='half_cauchy_alt.stan', seed=7878, refresh = 0)There are no warnings, and the sampling is as fast for the half-Cauchy nominal model.
samp<-as.array(fit_half_reparam)
res<-monitornew(samp[,,101:150], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
p1 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`sRhat`)) + geom_point() + ggtitle('Classic split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1)
p2 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`zsRhat`)) + geom_point() + ggtitle('Rank norm. split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1)
p3 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`fzsRhat`)) + geom_point() + ggtitle('Fold. rank split-Rhat') + geom_hline(yintercept = 1.005, linetype = 'dashed') + geom_hline(yintercept = 1)
grid.arrange(p1, p2, p3, nrow = 1)
All Rhats are below 1.01. Classic Rhat’s look also good even if the classic is not well defined for half-Cauchy distribution.
res<-monitornew(samp[,,101:150], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
p1 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`reff`)) + geom_point() + ggtitle('Classic relative efficiency (reff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + ylim(0,1)
p2 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`zsreff`)) + geom_point() + ggtitle('New relative efficiency (zsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + ylim(0,1)
p3 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=`fzsreff`)) + geom_point() + ggtitle('New relative efficiency (zsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + ylim(0,1)
p4 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=medsreff)) + geom_point() + ggtitle('Rel.eff. of median (medsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(limits = c(0,1))
p5 <- ggplot(data=as.data.frame(res), aes(x=1:50, y=madsreff)) + geom_point() + ggtitle('Rel.eff. of mad (madsreff)') + geom_hline(yintercept = c(0,1)) + geom_hline(yintercept = 0.1, linetype = 'dashed') + scale_y_continuous(limits = c(0,1))
grid.arrange(p1, p2, p3, blank, p4, p5, nrow = 2)
Result: Rank normalized r_eff’s have less variation than classic one which is not well defined for Cauchy. Based on rank normalized r_eff the alternative parameterization has much lower efficiency than the nominal parameterization with constraint <lower=0>.
We check the indicator function = and lp__
res<-monitornew(samp[,,151:152], warmup=0, digits_summary = 2, probs=0.5, print = FALSE);
cat('I: r_eff = ', round(res['I','zsreff'],2), '\n')I: r_eff = 0.78 cat('lp__: r_eff = ', round(res['lp__','zsreff'],2), '\n')lp__: r_eff = 0.24 We examine the relative efficiency in different parts of the posterior, by computing the relative efficiency for \(p(\hat{Q}_\alpha \le x < \hat{Q}_{\alpha+0.05})\), where \(\hat{Q}_\alpha\) is \(\alpha\)-quantile and \(\alpha \in (0,0.05,\ldots,0.95)\).
samp <- as.array(fit_half_reparam)
plot_local_reff(samp = samp, parami = 101, nalpha = 20)
We examine the relative efficiency in different parts of the posterior, by computing the relative efficiency for \(p(x < \hat{Q}_{\alpha+0.05})\), where \(\hat{Q}_\alpha\) is \(\alpha\)-quantile and \(\alpha \in (0,0.025,\ldots,0.975)\).
samp <- as.array(fit_half_reparam)
plot_quantile_reff(samp = samp, parami = 1, nalpha = 40)
The relative efficiency near zero is much worse than for the half-Cauchy with nominal parameterization and constraint <lower=0>.
Rank plot visualisation of x[15], which has the smallest relative efficiency 0.30 among x’s.
mcmc_hist(r_scale(samp[,,115]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/50)+0.5)
The rank plot is different from uniform, which is expected with lower relative efficiency, but would require a reference.
Rank plot visualisation of lp__, which has relative efficiency 0.24.
mcmc_hist(r_scale(samp[,,152]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/50)+0.5)
The rank plot is different from uniform, which is expected with lower relative efficiency, but would require a reference.
5 Hierarchical model: Eight Schools
The models are from Michael Betancourt’s case study on Diagnosing Biased Inference with Divergences.
5.1 A Centered Eight Schools model
writeLines(readLines("eight_schools_cp.stan"))data {
int<lower=0> J;
real y[J];
real<lower=0> sigma[J];
}
parameters {
real mu;
real<lower=0> tau;
real theta[J];
}
model {
mu ~ normal(0, 5);
tau ~ cauchy(0, 5);
theta ~ normal(mu, tau);
y ~ normal(theta, sigma);
}5.1.1 Centered parameterization default MCMC options
Run centered parameterization model with default MCMC options.
input_data <- read_rdump("eight_schools.data.R")
fit_cp <- stan(file='eight_schools_cp.stan', data=input_data,
iter=2000, chains=4, seed=483892929, refresh=0)Warning: There were 50 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupWarning: There were 2 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-lowWarning: Examine the pairs() plot to diagnose sampling problemsWe do observe divergent transitions, which is sensitive diagnostic for this model and we observer also low BFMI, but let’s assume we were investigating a MCMC algorithm without such diagnostic.
res<-monitornew(fit_cp, warmup=0, digits_summary = 2, probs=0.5, print = FALSE)
round(res[,c('sRhat','zsRhat','fzsRhat','reff','zsreff','fzsreff','medsreff','madsreff')],2) sRhat zsRhat fzsRhat reff zsreff fzsreff medsreff madsreff
mu 1.01 1.01 1.01 0.16 0.16 0.20 0.20 0.25
tau 1.02 1.02 1.00 0.07 0.05 0.22 0.09 0.24
theta[1] 1.00 1.00 1.01 0.28 0.26 0.16 0.25 0.19
theta[2] 1.00 1.00 1.01 0.27 0.26 0.21 0.26 0.21
theta[3] 1.01 1.01 1.01 0.21 0.21 0.22 0.24 0.25
theta[4] 1.00 1.00 1.02 0.29 0.28 0.29 0.27 0.23
theta[5] 1.01 1.01 1.01 0.22 0.22 0.29 0.23 0.33
theta[6] 1.00 1.00 1.01 0.28 0.27 0.23 0.24 0.24
theta[7] 1.00 1.00 1.01 0.23 0.22 0.20 0.25 0.22
theta[8] 1.00 1.00 1.01 0.31 0.28 0.19 0.26 0.30
lp__ 1.03 1.03 1.01 0.04 0.05 0.16 0.08 0.23Some rank-normalized split-Rhats are larger than 1.01. Relative efficiencies for tau and lp__ are less than 10%, which is worry-some and longer chains should be run.
We examine the relative efficiency of probabilities between quantiles (0, 0.05, …, 1) for tau.
samp <- as.array(fit_cp)
plot_local_reff(samp = samp, parami = 2, nalpha = 20)
We see that MCMC has difficulties in exploring small tau values. As the efficiency for estimating small tau values is practically zero, we may assume that we may miss substantial amount of posterior mass and get biased estimates.
We examine the relative efficiency of cumulative probabilities below quantiles (0.025, …, 0.975) for tau.
samp <- as.array(fit_cp)
plot_quantile_reff(samp = samp, parami = 2, nalpha = 40)
Rank plot visualisation of tau, which has the smallest relative efficiency 0.05 among parameters.
mcmc_hist(r_scale(samp[,,2]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/50)+0.5)
We observe clear problems.
Rank plot visualisation of lp__, which has relative efficiency 0.05.
mcmc_hist(r_scale(samp[,,11]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/50)+0.5)
We observe clear problems.
5.1.2 Centered parameterization with longer chains
Low r_eff can be sometimes compensated with longer chains. Let’s check 10 times longer chain.
fit_cp2 <- stan(file='eight_schools_cp.stan', data=input_data,
iter=20000, chains=4, seed=483892929, refresh=0)Warning: There were 1179 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupWarning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-lowWarning: Examine the pairs() plot to diagnose sampling problemsres<-monitornew(fit_cp2, warmup=0, digits_summary = 2, probs=0.5, print = FALSE)
round(res[,c('sRhat','zsRhat','fzsRhat','reff','zsreff','fzsreff','medsreff','madsreff')],2) sRhat zsRhat fzsRhat reff zsreff fzsreff medsreff madsreff
mu 1.00 1.00 1.00 0.07 0.07 0.17 0.03 0.06
tau 1.01 1.03 1.00 0.03 0.01 0.15 0.03 0.04
theta[1] 1.00 1.00 1.01 0.22 0.20 0.05 0.05 0.03
theta[2] 1.00 1.00 1.00 0.25 0.22 0.12 0.03 0.03
theta[3] 1.00 1.00 1.00 0.13 0.11 0.12 0.03 0.03
theta[4] 1.00 1.00 1.00 0.26 0.22 0.06 0.04 0.03
theta[5] 1.00 1.00 1.00 0.11 0.09 0.16 0.03 0.06
theta[6] 1.00 1.00 1.00 0.16 0.15 0.13 0.03 0.03
theta[7] 1.00 1.00 1.00 0.18 0.16 0.08 0.05 0.03
theta[8] 1.00 1.00 1.00 0.26 0.21 0.07 0.03 0.03
lp__ 1.02 1.03 1.01 0.01 0.01 0.01 0.04 0.05Some rank-normalized split-Rhats are still larger than 1.01. Relative efficiencies for tau and lp__ are around 1%. Drop in the relative efficiency when the the sample size was increased indicates serious problems in mixing.
We examine the relative efficiency of probabilities between quantiles (0, 0.05, …, 1) for tau.
samp <- as.array(fit_cp2)
plot_local_reff(samp = samp, parami = 2, nalpha = 50)
We see that MCMC has difficulties in exploring small tau values. As the efficiency for estimating small tau values is practically zero, we may assume that we may miss substantial amount of posterior mass and get biased estimates.
We examine the relative efficiency of cumulative probabilities below quantiles (0.01, …, 0.99) for tau.
samp <- as.array(fit_cp2)
plot_quantile_reff(samp = samp, parami = 2, nalpha = 100)
We examine also other parameters, like mu in the following plot.
samp <- as.array(fit_cp2)
plot_local_reff(samp = samp, parami = 1, nalpha = 50)
There are gaps of poor efficiency which indicates sticking of MCMC, but the sticking doesn’t occur for any specific range of values of mu.
We examine also other parameters, like mu in the following plot.
samp <- as.array(fit_cp2)
plot_quantile_reff(samp = samp, parami = 1, nalpha = 100)
Rank plot visualisation of tau, which has the smallest relative efficiency 0.01 among parameters.
mcmc_hist(r_scale(samp[,,2]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/100)+0.5)
We observe clear problems to sample small values of tau.
Rank plot visualisation of lp__, which has relative efficiency 0.01.
mcmc_hist(r_scale(samp[,,11]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/100)+0.5)
We observe clear problems sampling different energy levels.
5.1.3 Centered parameterization very long chains
And for further evidence let’s check 100 times longer chains than the default.
fit_cp3 <- stan(file='eight_schools_cp.stan', data=input_data,
iter=200000, chains=4, seed=483892929, refresh=0)Warning: There were 16912 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupWarning: There were 4 chains where the estimated Bayesian Fraction of Missing Information was low. See
http://mc-stan.org/misc/warnings.html#bfmi-lowWarning: Examine the pairs() plot to diagnose sampling problemsres<-monitornew(fit_cp3, warmup=0, digits_summary = 2, probs=0.5, print = FALSE)
round(res[,c('sRhat','zsRhat','fzsRhat','reff','zsreff','fzsreff','medsreff','madsreff')],2) sRhat zsRhat fzsRhat reff zsreff fzsreff medsreff madsreff
mu 1.00 1.00 1.00 0.04 0.03 0.01 0.01 0.01
tau 1.00 1.01 1.00 0.01 0.00 0.02 0.01 0.01
theta[1] 1.00 1.00 1.00 0.03 0.03 0.06 0.01 0.01
theta[2] 1.00 1.00 1.00 0.05 0.04 0.01 0.01 0.01
theta[3] 1.00 1.00 1.00 0.15 0.09 0.01 0.01 0.01
theta[4] 1.00 1.00 1.00 0.08 0.06 0.01 0.01 0.01
theta[5] 1.00 1.00 1.00 0.11 0.09 0.01 0.01 0.01
theta[6] 1.00 1.00 1.00 0.13 0.09 0.00 0.01 0.01
theta[7] 1.00 1.00 1.00 0.02 0.02 0.10 0.01 0.01
theta[8] 1.00 1.00 1.00 0.08 0.05 0.01 0.01 0.01
lp__ 1.01 1.01 1.01 0.00 0.00 0.00 0.01 0.01Some rank-normalized split-Rhats are still larger than 1.01. Relative efficiencies for tau and lp__ are around 0% and relative efficiencies for other parameters are also getting smaller. Drop in the relative efficiency when the the sample size was increased indicates serious problems in mixing.
We examine the relative efficiency of probabilities between quantiles (0, 0.05, …, 1) for tau.
samp <- as.array(fit_cp3)
plot_local_reff(samp = samp, parami = 2, nalpha = 100)
We see that MCMC has difficulties in exploring small tau values. As the efficiency for estimating small tau values is practically zero, we may assume that we may miss substantial amount of posterior mass and get biased estimates. It is good to note that the low efficiency for small tau values is due to too large step size, and the low efficiency for large tau values is partially due to too small step size, but even more importantly due to the random walk in energy space.
We examine the relative efficiency of cumulative probabilities below quantiles (0.01, …, 0.99) for tau.
samp <- as.array(fit_cp3)
plot_quantile_reff(samp = samp, parami = 2, nalpha = 100)
Rank plot visualisation of tau, which has the smallest relative efficiency 0.00 among parameters.
mcmc_hist(r_scale(samp[,,2]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/100)+0.5)
We observe clear problems to sample small values of tau and even with 100000 draws per chain, the plots don’t get crowded as traceplots would.
Rank plot visualisation of lp__, which has relative efficiency 0.01.
mcmc_hist(r_scale(samp[,,11]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/100)+0.5)
We observe clear problems sampling different energy levels, and the low efficiency for lp__ comes just from slow random walk (autocorrelation for chains 1 and 2 goes to zero after 1000 lags, and for chains 3 and 4 after 5000 lags).
5.2 Non-centered Eight Schools model
For hierarchical models, non-centered parameterization often works better
writeLines(readLines("eight_schools_ncp.stan"))data {
int<lower=0> J;
real y[J];
real<lower=0> sigma[J];
}
parameters {
real mu;
real<lower=0> tau;
real theta_tilde[J];
}
transformed parameters {
real theta[J];
for (j in 1:J)
theta[j] = mu + tau * theta_tilde[j];
}
model {
mu ~ normal(0, 5);
tau ~ cauchy(0, 5);
theta_tilde ~ normal(0, 1);
y ~ normal(theta, sigma);
}5.2.1 Non-centered parameterization default MCMC options
Run non-centered parameterization model with default options.
input_data <- read_rdump("eight_schools.data.R")
fit_ncp <- stan(file='eight_schools_ncp.stan', data=input_data,
iter=2000, chains=4, seed=483892929, refresh=0)Warning: There were 3 divergent transitions after warmup. Increasing adapt_delta above 0.8 may help. See
http://mc-stan.org/misc/warnings.html#divergent-transitions-after-warmupWarning: Examine the pairs() plot to diagnose sampling problemsWe still observe some divergent transitions with the default adapt_delta and no divergences with adapt_delta=0.95. Let’s analyze both samples.
res<-monitornew(fit_ncp, warmup=0, digits_summary = 2, probs=0.5, print = FALSE)
round(res[,c('sRhat','zsRhat','fzsRhat','reff','zsreff','fzsreff','medsreff','madsreff')],2) sRhat zsRhat fzsRhat reff zsreff fzsreff medsreff madsreff
mu 1 1 1 1.06 1.07 0.51 1.12 0.56
tau 1 1 1 0.66 0.57 0.70 0.83 0.77
theta_tilde[1] 1 1 1 1.27 1.29 0.49 1.25 0.53
theta_tilde[2] 1 1 1 1.33 1.35 0.54 1.23 0.58
theta_tilde[3] 1 1 1 1.33 1.34 0.47 1.28 0.58
theta_tilde[4] 1 1 1 1.29 1.32 0.50 1.26 0.63
theta_tilde[5] 1 1 1 1.25 1.25 0.40 1.11 0.53
theta_tilde[6] 1 1 1 1.28 1.30 0.47 1.27 0.60
theta_tilde[7] 1 1 1 1.13 1.15 0.56 1.22 0.59
theta_tilde[8] 1 1 1 1.48 1.50 0.52 1.41 0.58
theta[1] 1 1 1 0.95 1.02 0.60 1.05 0.69
theta[2] 1 1 1 1.07 1.11 0.59 1.14 0.67
theta[3] 1 1 1 1.05 1.13 0.63 1.19 0.63
theta[4] 1 1 1 1.12 1.16 0.61 1.12 0.66
theta[5] 1 1 1 1.15 1.19 0.55 1.25 0.63
theta[6] 1 1 1 0.99 1.05 0.65 1.11 0.69
theta[7] 1 1 1 0.97 1.07 0.62 1.12 0.76
theta[8] 1 1 1 1.12 1.18 0.64 1.29 0.73
lp__ 1 1 1 0.36 0.38 0.64 0.47 0.84All Rhats are close to 1, and relative efficiencies are good.
We examine the relative efficiency of probabilities between quantiles (0, 0.05, …, 1) for tau.
samp <- as.array(fit_ncp)
plot_local_reff(samp = samp, parami = 2, nalpha = 20)
Small tau values are still more difficult to explore, but the relative efficiency is in a good range.
We examine the relative efficiency of probabilities below quantiles (0.025, …, 0.975) for tau.
samp <- as.array(fit_ncp)
plot_quantile_reff(samp = samp, parami = 2, nalpha = 40)
Rank plot visualisation of tau, which has the smallest relative efficiency 0.57 among parameters.
mcmc_hist(r_scale(samp[,,2]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/50)+0.5)
We can see problems in sampling small values of tau.
Rank plot visualisation of lp__, which has relative efficiency 0.01.
mcmc_hist(r_scale(samp[,,19]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/50)+0.5)
We observe clear problems sampling different energy levels.
5.2.2 Non-centered parameterization default MCMC options plus adapt_delta=0.95
Next we examine the same model but with higher adapt_delta=0.95.
fit_ncp2 <- stan(file='eight_schools_ncp.stan', data=input_data,
iter=2000, chains=4, control=list(adapt_delta=0.95), seed=483892929, refresh=0)We get zero divergences with adapt_delta=0.95.
res<-monitornew(fit_ncp2, warmup=0, digits_summary = 2, probs=0.5, print = FALSE)
round(res[,c('sRhat','zsRhat','fzsRhat','reff','zsreff','fzsreff','medsreff','madsreff')],2) sRhat zsRhat fzsRhat reff zsreff fzsreff medsreff madsreff
mu 1 1 1 1.15 1.17 0.53 1.13 0.61
tau 1 1 1 0.83 0.62 0.77 0.82 0.70
theta_tilde[1] 1 1 1 1.09 1.11 0.54 1.14 0.63
theta_tilde[2] 1 1 1 1.16 1.17 0.45 1.05 0.56
theta_tilde[3] 1 1 1 1.23 1.23 0.45 1.16 0.56
theta_tilde[4] 1 1 1 1.26 1.28 0.45 1.26 0.58
theta_tilde[5] 1 1 1 1.15 1.16 0.49 1.17 0.56
theta_tilde[6] 1 1 1 1.29 1.29 0.52 1.19 0.59
theta_tilde[7] 1 1 1 1.10 1.11 0.55 1.09 0.64
theta_tilde[8] 1 1 1 1.19 1.19 0.53 1.22 0.60
theta[1] 1 1 1 0.99 1.07 0.66 1.13 0.72
theta[2] 1 1 1 1.24 1.29 0.58 1.21 0.66
theta[3] 1 1 1 1.07 1.13 0.68 1.07 0.73
theta[4] 1 1 1 1.24 1.27 0.67 1.27 0.73
theta[5] 1 1 1 1.17 1.20 0.72 1.08 0.77
theta[6] 1 1 1 1.11 1.18 0.63 1.21 0.64
theta[7] 1 1 1 1.18 1.22 0.68 1.17 0.74
theta[8] 1 1 1 1.03 1.10 0.69 1.15 0.71
lp__ 1 1 1 0.44 0.44 0.64 0.51 0.70All Rhats are close to 1, and relative efficiencies are good and slightly better than with the default adapt_delta.
We examine the relative efficiency of probabilities between quantiles (0, 0.05, …, 1) for tau.
samp <- as.array(fit_ncp2)
plot_local_reff(samp = samp, parami = 2, nalpha = 20)
Small tau values are still more difficult to explore, but the relative efficiency is in a good range.
Check with a finer resolution
samp <- as.array(fit_ncp2)
plot_local_reff(samp = samp, parami = 2, nalpha = 50)
Also with a finer resolution the efficiency of MCMC is good for different values of tau.
We examine the relative efficiency of probabilities below quantiles (0.025, …, 0.975) for tau.
samp <- as.array(fit_ncp2)
plot_quantile_reff(samp = samp, parami = 2, nalpha = 40)
Rank plot visualisation of tau, which has the smallest relative efficiency 0.62 among parameters.
mcmc_hist(r_scale(samp[,,2]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/100)+0.5)
Higher adapt_delta seems to have improved sampling of small values.
Rank plot visualisation of lp__, which has relative efficiency 0.44.
mcmc_hist(r_scale(samp[,,19]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/100)+0.5)
We observe clear problems sampling different energy levels.
5.2.3 Non-centered parameterization default MCMC options plus adapt_delta=0.95 and longer chains
If in doubt, we can run longer chains.
input_data <- read_rdump("eight_schools.data.R")
start_time <- Sys.time();
fit_ncp3 <- stan(file='eight_schools_ncp.stan', data=input_data,
iter=20000, chains=4, control=list(adapt_delta=0.954), seed=483892929, refresh=0)
end_time <- Sys.time();end_time - start_timeres<-monitornew(fit_ncp3, warmup=0, digits_summary = 2, probs=0.5, print = FALSE)
round(res[,c('sRhat','zsRhat','fzsRhat','reff','zsreff','fzsreff','medsreff','madsreff')],2) sRhat zsRhat fzsRhat reff zsreff fzsreff medsreff madsreff
mu 1 1 1 1.09 1.10 0.51 1.11 0.62
tau 1 1 1 0.87 0.65 0.78 0.86 0.75
theta_tilde[1] 1 1 1 1.21 1.21 0.54 1.15 0.63
theta_tilde[2] 1 1 1 1.22 1.23 0.49 1.19 0.56
theta_tilde[3] 1 1 1 1.20 1.20 0.48 1.17 0.55
theta_tilde[4] 1 1 1 1.23 1.23 0.50 1.23 0.58
theta_tilde[5] 1 1 1 1.25 1.26 0.53 1.18 0.62
theta_tilde[6] 1 1 1 1.23 1.23 0.48 1.21 0.56
theta_tilde[7] 1 1 1 1.15 1.16 0.53 1.12 0.62
theta_tilde[8] 1 1 1 1.31 1.31 0.47 1.22 0.56
theta[1] 1 1 1 1.07 1.12 0.66 1.10 0.71
theta[2] 1 1 1 1.19 1.21 0.65 1.15 0.70
theta[3] 1 1 1 1.04 1.08 0.63 1.11 0.70
theta[4] 1 1 1 1.13 1.15 0.64 1.14 0.68
theta[5] 1 1 1 1.21 1.21 0.65 1.14 0.72
theta[6] 1 1 1 1.15 1.19 0.64 1.13 0.70
theta[7] 1 1 1 1.12 1.12 0.68 1.08 0.74
theta[8] 1 1 1 1.07 1.11 0.60 1.08 0.69
lp__ 1 1 1 0.41 0.41 0.64 0.51 0.75All Rhats are close to 1, and relative efficiencies are good and slightly better than with the default adapt_delta. Relative efficiency is similar as for shorter chain, which means that running longer chains is giving us higher effective sample size.
We examine the relative efficiency of probabilities between quantiles (0, 0.05, …, 1) for tau.
samp <- as.array(fit_ncp3)
plot_local_reff(samp = samp, parami = 2, nalpha = 100)
Small tau values are still more difficult to explore, but the relative efficiency is in a good range.
We examine the relative efficiency of probabilities below quantiles (0.01, …, 0.99) for tau.
samp <- as.array(fit_ncp2)
plot_quantile_reff(samp = samp, parami = 2, nalpha = 100)
Rank plot visualisation of tau, which has the smallest relative efficiency 0.64 among parameters.
mcmc_hist(r_scale(samp[,,2]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/200)+0.5)
Higher adapt_delta seems to have improved sampling of small values.
Rank plot visualisation of lp__, which has relative efficiency 0.41.
mcmc_hist(r_scale(samp[,,19]), breaks = seq(0,prod(dim(samp)[1:2]),by=prod(dim(samp)[1:2])/200)+0.5)
We observe clear problems sampling different energy levels.
6 NUTS and effective sample size
We have already seen that the relative efficiency of HMC/NUTS can be higher than with independent draws. The next example illustrates interesting relative efficiency phenomena due to properties of HMC/NUTS algorithm.
writeLines(readLines("normal.stan"))data {
int<lower=1> J;
}
parameters {
vector[J] x;
}
model {
x ~ normal(0, 1);
}fit_n <- stan(file='normal.stan', data=data.frame(J=16),
iter=20000, chains=4, seed=483892929, refresh=0)samp<-as.array(fit_n)
res<-monitornew(samp, warmup=0, digits_summary = 2, probs=0.5, print = FALSE)
round(res[,c('sRhat','zsRhat','fzsRhat','reff','zsreff','fzsreff','medsreff','madsreff')],2) sRhat zsRhat fzsRhat reff zsreff fzsreff medsreff madsreff
x[1] 1 1 1 2.53 2.55 0.39 2.01 0.45
x[2] 1 1 1 2.67 2.68 0.41 2.01 0.48
x[3] 1 1 1 2.54 2.56 0.41 2.05 0.46
x[4] 1 1 1 2.53 2.54 0.40 2.10 0.48
x[5] 1 1 1 2.49 2.50 0.41 2.00 0.47
x[6] 1 1 1 2.51 2.53 0.42 2.09 0.50
x[7] 1 1 1 2.49 2.50 0.39 2.20 0.45
x[8] 1 1 1 2.64 2.65 0.40 2.01 0.45
x[9] 1 1 1 2.55 2.56 0.41 2.04 0.49
x[10] 1 1 1 2.62 2.63 0.40 2.03 0.47
x[11] 1 1 1 2.54 2.55 0.41 1.99 0.48
x[12] 1 1 1 2.50 2.52 0.40 2.01 0.47
x[13] 1 1 1 2.59 2.60 0.41 2.08 0.49
x[14] 1 1 1 2.49 2.51 0.41 2.11 0.48
x[15] 1 1 1 2.58 2.59 0.40 2.13 0.46
x[16] 1 1 1 2.68 2.69 0.41 2.06 0.48
lp__ 1 1 1 0.36 0.35 0.53 0.42 0.59The relative efficiency for all \(x\) is larger than 2. The relative efficiency for lp__ is however only 0.35. If look at the all Stan examples in this notebook, we see that the relative efficiency for lp__ is always below 0.5. lp__ correlates strongly with the total energy in HMC, and that total energy is sampled using random walk proposal once per iteration and thus it’s likely that lp__ has some random walk behavior leading to autocorrelation and relative efficiency below 1. On the other hand, no-U-turn behavior and sampling from the trajectory makes the Markov chain to be antithetic with negative odd lag correlation and relative efficiency higher than 1.
Let’s check the local relative efficiency
plot_local_reff(samp = samp, parami = 1, nalpha = 100)
The relative efficiency for probability estimate for a small interval is close to 1 with a slight drop in the tails. This is a good result, but far from the relative efficiency for mean estimate.
Let’s check the relative efficiency of cumulative probabilities
plot_quantile_reff(samp = samp, parami = 1, nalpha = 100)
The total energy of HMC should affect how far in the tails a chain in one iteration can go. Far tails of the target have high energy, and only chains with high total energy can then reach there. This will suggest that, the random walk in total energy would cause random walk in variance of \(x\). Let’s check the second moment of \(x\).
samp<-as.array(fit_n)^2
res<-monitornew(samp, warmup=0, digits_summary = 2, probs=0.5, print = FALSE)
round(res[,c('sRhat','zsRhat','fzsRhat','reff','zsreff','fzsreff','medsreff','madsreff')],2) sRhat zsRhat fzsRhat reff zsreff fzsreff medsreff madsreff
x[1] 1 1 1 0.36 0.39 0.47 0.45 0.60
x[2] 1 1 1 0.37 0.41 0.50 0.48 0.63
x[3] 1 1 1 0.38 0.41 0.48 0.46 0.60
x[4] 1 1 1 0.37 0.40 0.47 0.48 0.60
x[5] 1 1 1 0.37 0.41 0.47 0.47 0.59
x[6] 1 1 1 0.37 0.42 0.47 0.50 0.61
x[7] 1 1 1 0.35 0.39 0.46 0.46 0.58
x[8] 1 1 1 0.36 0.40 0.46 0.46 0.57
x[9] 1 1 1 0.38 0.41 0.48 0.48 0.60
x[10] 1 1 1 0.36 0.40 0.47 0.47 0.61
x[11] 1 1 1 0.37 0.41 0.48 0.48 0.59
x[12] 1 1 1 0.36 0.40 0.47 0.47 0.57
x[13] 1 1 1 0.36 0.41 0.48 0.49 0.62
x[14] 1 1 1 0.37 0.41 0.47 0.48 0.61
x[15] 1 1 1 0.37 0.40 0.47 0.46 0.60
x[16] 1 1 1 0.36 0.41 0.47 0.48 0.57
lp__ 1 1 1 0.39 0.35 0.53 0.42 0.58Mean of the relative efficiencies for \(x[j]^2\) is 0.4, which is quite close to the relative efficiency for lp__. This is not that surprising as the potential energy in normal model is relative to \(\sum_{j=1}^J x_j^2\).
Let’s check the local relative efficiency
plot_local_reff(samp = samp, parami = 1, nalpha = 100)
The relative efficiency is mostly a bit below 1, but for the right tail of \(x[1]^2\) the relative efficiency drops. This likely due to only some iterations have high enough total energy to obtain draws from the high energy part of the tail.
Let’s check the relative efficiency of cumulative probabilities
plot_quantile_reff(samp = samp, parami = 1, nalpha = 100)
We can see the correlation between lp__ and magnitude of \(x[1]\) in the following plot.
samp<-as.array(fit_n)
qplot(as.vector(samp[,,17]),abs(as.vector(samp[,,1]))) + xlab('lp__') + ylab('x[1]')
Low lp__ corresponds to high energy and more variation in \(x\), and high lp__ corresponds to low energy and small variation in \(x\).
Finally \(\sum_{j=1}^J x_j^2\) is perfectly correlated with lp__.
qplot(as.vector(samp[,,17]),as.vector(apply(samp[,,1:16]^2,c(1,2),sum))) + xlab('lp__') + ylab('sum(x^2)')
This shows that even if we get high relative efficiency estimates for parameters, it is important to look at the relative efficiency of lp__ as it can indicate problems of sampling in tails, which would affect the accuracy of variance and tail quantile estimates.
7 Local efficiency for distributional approximations
We note that we can plot similar local relative efficiency plots also for distributional approximations (like Laplace and ADVI approximations), using Pareto smoothed importance sampling based effective sample size estimates (Vehtari et al, 2017). There will be another case study for this.
8 Original Computing Environment
writeLines(readLines(file.path(Sys.getenv("HOME"), ".R/Makevars")))CXXFLAGS=-O3 -mtune=native -march=native -Wno-unused-variable -Wno-unused-function
CXXFLAGS+=-flto -ffat-lto-objects -Wno-unused-local-typedefs
CXXFLAGS+=-std=c++11
CFLAGS+=-O3devtools::session_info("rstan")Session info ------------------------------------------------------------- setting value
version R version 3.5.1 (2018-07-02)
system x86_64, linux-gnu
ui X11
language en_GB:en
collate en_US.UTF-8
tz Europe/Helsinki
date 2018-08-16 Packages ----------------------------------------------------------------- package * version date source
assertthat 0.2.0 2017-04-11 CRAN (R 3.5.1)
BH 1.66.0-1 2018-02-13 CRAN (R 3.5.1)
cli 1.0.0 2017-11-05 CRAN (R 3.5.1)
colorspace 1.3-2 2016-12-14 CRAN (R 3.5.1)
crayon 1.3.4 2017-09-16 CRAN (R 3.5.1)
dichromat 2.0-0 2013-01-24 CRAN (R 3.5.1)
digest 0.6.15 2018-01-28 CRAN (R 3.5.1)
fansi 0.2.3 2018-05-06 CRAN (R 3.5.1)
ggplot2 * 3.0.0 2018-07-03 CRAN (R 3.5.1)
glue 1.3.0 2018-07-17 CRAN (R 3.5.1)
graphics * 3.5.1 2018-07-03 local
grDevices * 3.5.1 2018-07-03 local
grid * 3.5.1 2018-07-03 local
gridExtra * 2.3 2017-09-09 CRAN (R 3.5.1)
gtable 0.2.0 2016-02-26 CRAN (R 3.5.1)
inline 0.3.15 2018-05-18 CRAN (R 3.5.1)
labeling 0.3 2014-08-23 CRAN (R 3.5.1)
lattice 0.20-35 2017-03-25 CRAN (R 3.5.1)
lazyeval 0.2.1 2017-10-29 CRAN (R 3.5.1)
magrittr 1.5 2014-11-22 CRAN (R 3.5.1)
MASS 7.3-50 2018-04-30 CRAN (R 3.5.0)
Matrix 1.2-14 2018-04-09 CRAN (R 3.5.1)
methods * 3.5.1 2018-07-03 local
mgcv 1.8-24 2018-06-18 CRAN (R 3.5.0)
munsell 0.5.0 2018-06-12 CRAN (R 3.5.1)
nlme 3.1-137 2018-04-07 CRAN (R 3.5.1)
pillar 1.3.0 2018-07-14 CRAN (R 3.5.1)
plyr 1.8.4 2016-06-08 CRAN (R 3.5.1)
R6 2.2.2 2017-06-17 CRAN (R 3.5.1)
RColorBrewer 1.1-2 2014-12-07 CRAN (R 3.5.1)
Rcpp 0.12.18 2018-07-23 CRAN (R 3.5.1)
RcppEigen 0.3.3.4.0 2018-02-07 CRAN (R 3.5.1)
reshape2 1.4.3 2017-12-11 CRAN (R 3.5.1)
rlang 0.2.1 2018-05-30 CRAN (R 3.5.1)
rstan * 2.17.3 2018-01-20 CRAN (R 3.5.1)
scales 0.5.0 2017-08-24 CRAN (R 3.5.1)
StanHeaders * 2.17.2 2018-01-20 CRAN (R 3.5.1)
stats * 3.5.1 2018-07-03 local
stats4 3.5.1 2018-07-03 local
stringi 1.2.4 2018-07-20 CRAN (R 3.5.1)
stringr 1.3.1 2018-05-10 CRAN (R 3.5.1)
tibble 1.4.2 2018-01-22 CRAN (R 3.5.1)
tools 3.5.1 2018-07-03 local
utf8 1.1.4 2018-05-24 CRAN (R 3.5.1)
utils * 3.5.1 2018-07-03 local
viridisLite 0.3.0 2018-02-01 CRAN (R 3.5.1)
withr 2.1.2 2018-03-15 CRAN (R 3.5.1)